function F = solve_eqm_outer(param, const)
%% Initial Equilibrium

    % Initial n(q)
    q_indicator = sparse(param.Nf, param.Q);
    q_indicator(:,1) = 1; %start from the lowest quality  
    nq_new = cal_nq(const, q_indicator); 

    % Initial integrals
    integ = cal_integral(const, q_indicator, ones(param.Nf,param.Q), ones(param.Nf,param.Q), ones(param.Nf,param.Q));   

    % Initial Pi(q,E,I)
    eqm.Pi00 = ones(1,param.Q);
    eqm.Pi01 = ones(1,param.Q);
    eqm.Pi10 = ones(1,param.Q);
    eqm.Pi11 = ones(1,param.Q);
    
    % Settings
    tol = 1e-12;
    maxIter = 100;
    maxHistory = 5;
    maxRepeat = maxHistory+1;
    
    % Placeholder/initialize (check repeat)
    q_opt_index_history = zeros(param.Nf,maxHistory); 
    switch_id_history = zeros(param.Nf,maxHistory+1);
    switch_qindex_history = zeros(param.Nf,maxHistory+1);     
    fix_flag = 0;
    repeat = zeros(1,maxHistory);
    
    % Figure
    figure; 
    
    % Iteration
    err = 1;
    iter = 0;
    while (err >= tol && iter < maxIter)

        % save previous H quality choice
        if iter > 1
            q_opt_index_history = [q_opt_index, q_opt_index_history(:,[1:end-1])];
        end   
        
        % save previous H switch id/qindex
        if max(repeat(2:end)) > 0           
            temp_switch_id = zeros(param.Nf,1);
            temp_switch_id(switch_id) = switch_id;
            switch_id_history = [temp_switch_id, switch_id_history(:,[1:end-1])];            
            temp_switch_qindex = zeros(param.Nf,1);
            temp_switch_qindex(switch_id) = switch_qindex;
            switch_qindex_history = [temp_switch_qindex, switch_qindex_history(:,[1:end-1])];
        end
        
        % new guess
        Pi00 = eqm.Pi00;
        Pi01 = eqm.Pi01;
        Pi10 = eqm.Pi10;
        Pi11 = eqm.Pi11;
        nq = nq_new;
        iter = iter + 1;
                
        % solve equilibrium
        eqm = solve_eqm_inner(param, const, integ); 
        
        % choose new quality (grid search)
        [q_opt_index, ~] = grid_search(param, eqm.EPiQ0, eqm.EPiQ1, eqm.ExportProbQ, eqm.ExportCutoffQ);        

        % fix quality after maxRepeat
        if max(repeat(2:end)) >= maxRepeat 
            if fix_flag == 0
                % find id and qindex
                fix_id = nonzeros(max(switch_id_history,[],2));
                fix_qindex_history = switch_qindex_history(fix_id,:);
                fix_qindex = max(fix_qindex_history,[],2); 
                % set flag
                fix_flag = 1;
                repeat = repeat*0;
            elseif fix_flag > 0
                % find id and qindex
                fix_id2 = nonzeros(max(switch_id_history,[],2));
                fix_qindex_history2 = switch_qindex_history(fix_id2,:);
                fix_qindex2 = max(fix_qindex_history2,[],2);  
                % append
                fix_id = [fix_id; fix_id2];
                fix_qindex = [fix_qindex; fix_qindex2];
                fix_qindex_history = [fix_qindex_history; fix_qindex_history2];
                % set flag
                fix_flag = fix_flag + 1;
                repeat = repeat*0;
            end
        end
        if fix_flag > 0
            q_opt_index(fix_id) = fix_qindex;                
        end          
        
        % update qindicator
        q_indicator = cal_indicator(param, q_opt_index); 

        % update integrals
        integ = cal_integral(const, q_indicator, eqm.ExportProbQ, eqm.ImportProbQ0, eqm.ImportProbQ1); 
        
        % update n(q)
        nq_new = cal_nq(const, q_indicator);          
        
        % check convergence 
        err1 = max(abs(nq_new-nq));
        err2 = max(abs(eqm.Pi00-Pi00));
        err3 = max(abs(eqm.Pi01-Pi01));        
        err4 = max(abs(eqm.Pi10-Pi10));
        err5 = max(abs(eqm.Pi11-Pi11));
        err = max([err1, err2, err3, err4, err5]);
        
            % visualize progress
            scatter(const.Qgrid, nq_new, 'Marker', '.')
            xlabel('q')  
            title(strcat('iter=',num2str(iter)))
            drawnow limitrate
        
        % firm types changing quality
        switch_id = find(q_opt_index~=q_opt_index_history(:,1));
        switch_qindex = q_opt_index(switch_id);
        switch_num = size(switch_id,1);
        
        % repeat history
        for i = 1:maxHistory
            repeat(i) = repeat(i)*min(q_opt_index == q_opt_index_history(:,i)) + min(q_opt_index == q_opt_index_history(:,i));
        end                                  
        
            % print
            fprintf('iter%.0f: %.0f/%.0f types change quality', iter, switch_num, param.Nf); 
            if switch_num > 0
                if max(repeat(2:end)) > 0
                    fprintf(', repeat ='); fprintf(' %.0f', repeat); 
                end
                if switch_num <= 5
                    fprintf(', id ='); fprintf(' %.0f', switch_id);
                    fprintf(', qindex ='); fprintf(' %.0f', switch_qindex);
                end 
            elseif switch_num == 0
                fprintf(', err = %e', err);
            end
            fprintf('\n')

    end
    close
    

%% Store   
    
    eqm.iter_outer      = iter;
    eqm.err_outer       = err;
    if fix_flag > 0    
        eqm.fix_id               = fix_id;
        eqm.fix_qindex           = fix_qindex; 
        eqm.fix_qindex_history   = fix_qindex_history;
    else
        eqm.fix_id               = nan;
        eqm.fix_qindex           = nan;         
        eqm.fix_qindex_history   = nan;
    end  
    
    eqm.q_index         = q_opt_index;
    eqm.q_indicator     = q_indicator;
    eqm.integ           = integ;  
    
    eqm.nq               = nq_new;
    eqm.first_nonzero    = find(eqm.nq, 1, 'first');
    eqm.last_nonzero     = find(eqm.nq, 1, 'last');  
    eqm.first_nonzero_q  = const.Qgrid(eqm.first_nonzero);
    eqm.last_nonzero_q   = const.Qgrid(eqm.last_nonzero);  
    
    eqm.density_all     = const.density_omega;  
    eqm.density_exp     = const.density_omega.*full(sum(eqm.ExportProbQ.*eqm.q_indicator,2));
    eqm.density_nonexp  = const.density_omega - eqm.density_exp;
    
    
%% 
    % import probability of exporters
    eqm.prob_exp_imp        = full(sum(eqm.ImportProbQ1.*eqm.q_indicator, 2));
    eqm.prob_exp_nonimp     = full(sum((1-eqm.ImportProbQ1).*eqm.q_indicator, 2));
    
    % import probability of non-exporters
    eqm.prob_nonexp_imp     = full(sum(eqm.ImportProbQ0.*eqm.q_indicator, 2));
    eqm.prob_nonexp_nonimp  = full(sum((1-eqm.ImportProbQ0).*eqm.q_indicator, 2));
    
        % check
        check1 = eqm.prob_exp_imp + eqm.prob_exp_nonimp;
        check0 = eqm.prob_nonexp_imp + eqm.prob_nonexp_nonimp;
        
    % density of importer/non-importer
    eqm.density_imp =  eqm.density_exp.*eqm.prob_exp_imp + eqm.density_nonexp.*eqm.prob_nonexp_imp;
    eqm.density_nonimp = const.density_omega - eqm.density_imp;

    
%% Determine Quintile Group for each q

    % Quintileflag
    quintileflag = zeros(param.d, param.Q);
    quintilecutoff = zeros(1, param.d+1);
    nq_cum = cumsum(eqm.nq);
    for i = 1: param.d
        [~,quintilecutoff(i+1)] = min(abs(nq_cum - 0.2*i));
        quintileflag(i,quintilecutoff(i)+1:quintilecutoff(i+1)) = 1;            
    end
    quintileflag(1,(eqm.nq==0))=0;
    nd = sum(repmat(eqm.nq, param.d,1).*quintileflag,2)';  
        
    % Store
    eqm.nd                = nd;
    eqm.quintilecutoff    = quintilecutoff;
    eqm.quintileflag      = quintileflag;    
    
    
%% Indicators for Quintile Group for each firm type

    % Indicators for ex-ante quintile
    eqm.Qindicator = eqm.q_indicator*(eqm.quintileflag)'; %Nf by 5 matrix
    
    
%% Return

    F = eqm;
 
    
end